Impurities in a Biased Graphene Bilayer 
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We study the problem of impurities and mid-gap states in a biased graphene bilayer. We show 
that the properties of the bound states, such as localization lengths and binding energies, can be 
controlled externally by an electric field effect. Moreover, the band gap is renormalized and impurity 
bands are created at finite impurity concentrations. Using the coherent potential approximation we 
calculate the electronic density of states and its dependence on the applied bias voltage. 

PACS numbers: 81.05.Uw 73.21.Ac 73.20.Hb 



The discovery of single layer graphene and its 
further experimental characterization demonstrating un- 
conventional metallic properties 0, have attracted 
a lot of attention in the condensed matter community. 
The graphene bilayer, which is made out of two stacked 
graphene planes, is a particularly interesting form of two- 
dimensional (2D) carbon because of its unusual physics 
that has its origins on the peculiar band structure. At 
low energies and long wavelengths it can be described in 
terms of massive, chiral, Dirac particles. Previous stud- 
ies of the graphene bilayer have focused on the integer 
quantum Hall effect [1, 0] , the effect of electron-electron 
interactions @ , and transport properties @, H, B 03 U | ■ 
An important property of this system is that when an 
electric field is applied between the two graphene lay- 
ers, an electronic gap opens in the spectrum 0, IH, EH; 
the resulting system is called the biased graphene bilayer 
(BGB). In contrast, the single layer graphene does not 
have this unique property since to generate a gap in the 
electronic spectrum the sub-lattice symmetry of the hon- 
eycomb lattice has to broken, and this is a costly ener- 
getic process. 

Ohta et al. have studied graphene films on SiC sub- 
strates using angle-resolved photoemission spectroscopy 
(ARPES) [14|, and spectra reminiscent of the BGB were 
uncovered. Nevertheless, the graphene films are heavily 
doped by the SiC and the chemical potential and gap 
value cannot be controlled independently. More recently, 
Castro et al. reported the first observation of a tunable 
electronic gap in a BGB through magneto-transport mea- 
surements of micro-mechanically cleaved graphene on a 
Si02 substrate [15|- By chemically doping the upper 
graphene layer with NH3, it was shown that the BGB 
behaves as a semiconductor with a tunable electronic gap 
that can be changed from zero to a value as large as 0.3 
eV by using fields of <, 1 V/nm (below the electric break- 
down of Si02). The importance of controlling both the 
chemical potential and the value of the electronic gap in 
a semiconductor cannot be overstated since it can open 
doors for a large number of applications, from transistors 
[lit to photo-detectors and lasers tunable by the electric 
field effect. 

In order to understand and control the electronic prop- 



erties of the BGB it is important to understand the effects 
of the unavoidable disorder. In this letter we show that 
bound states exist for arbitrary weak impurity potentials, 
and that their properties, such as binding energies and 
localization lengths, can be externally controlled with a 
gate bias. Moreover, we obtain the wave-functions of 
the mid-gap states, from which we derive a simple crite- 
rion for when the overlap between wavefunctions becomes 
important. This overlap results in band gap renormal- 
ization and possibly band tails extending into the gap 
region, as in the case of ordinary heavily doped semi- 
conductors [l7j |. or impurity bands for deep impurities. 
Unlike ordinary semiconductors, the electronic density 
of states can be completely controlled via the electric 
field effect. The impurity interaction problem is studied 
within the coherent potential approximation (CPA). 

The Model. The low-energy effective bilayer Hamilto- 
nian has the form 12, IH (we use units such that % = 1): 



fa* -k + V/2 t x (l + a z )/2\ 
Wo(k) - {t ± (l + a z )/2 a-k-V/2 J' (1) 

where k = {k Xl k v ) is the 2D momentum measured rel- 
ative to the K-point in the Brillouin zone (BZ), V is 
the potential energy difference between the two planes, 
t± « 0.35 eV is the inter-layer hopping energy, and 
Oi (i = x, y, z) are Pauli matrices. We choose units 
such that the Fermi-Dirac velocity, up, is set to unity 
(vf = Ma/ 2 where t w 3eV is the nearest neighbor hop- 
ping energy, and a w 1.4 A is the lattice spacing). The 
corresponding spinor has weight on the different sublat- 
tices according to * = (^ A1 , ^ B1 , ip A2 , Vb2) T @- Solv- 
ing for the spectrum of (JTJ) one finds two pairs of electron- 
hole symmetric bands: 



k + V z /A + t z j2 ±J(V 2 + t\)V + 4/4. (2) 



For the two bands closest to zero energy near the band 
edge (valence and conduction bands) one can write the 
spectrum as: E±-(k) « ±[E s /2 + (k - fc g ) 2 /(2m g )]. 
Here E s = Vt±/ ' \JV 2 + 1\ is the energy gap, fc g = 
V/2y/l + t'^/CV 2 + tj_) is a momentum shift, and m g is 
an effective mass. Due to the non-zero k g the system is 
effectively one-dimensional (ID) near the band edge 
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Dirac delta potentials. True bound states must lie in- 
side of the gap so that their energies fulfill |e| < -Eg/2. 
For a single local impurity the Green's function can be 
written as G = G° + G°TG° using the standard t- 
matrix approach (l8j . where the bare Green's function 
is G° = (ui — 7io) _1 - The Fourier transform of a Dirac 
delta potential is U/N (N being the number of unit cells) 
leading to T = (U/N)/(l - U(f aj ). Here a = (A, B) and 
j = (1, 2) label the lattice site of the impurity. The quan- 
tity G aj (to) = J2k ^aj a 3 (w, k)/iV is the local propagator 
at the impurity site and the momentum sum is over the 
whole BZ. Bound states are then identified by the addi- 
tional poles of the Green's function due to the potential, 
these are given by G QJ (e) = 1/U. For energies inside of 
the gap we find: 

-o V/2 -u r ( A 4 



,m 4 + <y 2 /4- 
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where M 2 = yJV 2 t\/A - lu 2 (V 2 + t\), and A (w 7 eV) is 
a high energy cut-off Q ■ The corresponding expressions 
in plane 2 are obtained by the substitution V — > — V. 
From this we conclude that a Dirac delta potential al- 
ways generates a bound state since G diverges as the 
band edge is approached (where M — > 0). The depen- 
dence on the cut-off (except for the overall scale) is rather 
weak so that the linear in-plane approximation to the 
spectrum should be a good approximation as in the case 
of graphene For a given strength of the potential 

U, there are four different bound state energies depend- 
ing on which lattice site it is sitting on. In Fig. [1] we 
show the binding energy as a function of U and V for 
the deepest bound state. In the limit of U — > oo the 
particle-hole symmetry of the bound state energies is re- 
stored. We will only consider attractive potentials in this 
work, analogous results will hold for repulsive potentials 
because of the particle-hole symmetry of the model. For 
smaller values of the potential (\U\ <C A) the binding en- 
ergy measured from the band edge: E\> = E g /2 — e grows 
as U 2 and the states are weakly bound. For example, for 
V = 40mcV and U < leV one finds E h < 4 x lCT 4 £; g . 

Angular momentum states. For any potential with 
cylindrical symmetry it is useful to classify the eigen- 
states according to their angular momentum m. The 
real-space version of (JT|) in polar coordinates is ob- 
tained by the substitutions k x ± ik y — > —i{d x ± id y ) — > 
—ie ±zlp (dr ± id^/r). Adding a potential that only de- 
pends on the radial coordinate r one can (in analogy with 
the usual Dirac equation |20() construct an angular mo- 
mentum operator that commutes with the Hamiltonian. 
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FIG. 1: (Color online) Left: Bound state binding energies 
(in units of the gap E s ) for a Dirac delta potential of strength 
~U for different bias V. Right: Binding energy of a potential 
well of range R — 10 a and strength 7 = 71 = 72 (see text) 
for different angular momentum m. 



The angular (ip) dependence of the angular momentum 
m eigenstates are those of the vector: 



i Q , m (<y2) = e imip [l, e-' v e- ia7T/2 , 1, e ilp e ia7r ^ 2 ] 



(4) 



The parameter a is introduced for later convenience. If 
the potential generates bound states inside of the gap 
these states decay exponentially: ~ r 7 e _Kr . Assuming 
that the potential decays fast enough, the asymptotic 
behavior of ([1} implies that the allowed values of k are 



k± = vMe 2 + V 2 /A)±iM 2 , (5) 
so that weakly bound states have a localization length 



l~[2k s /(Vt ± )]^/E s /E h , 



(6) 



that diverges as the band edge is approached and de- 
creases with increasing bias voltage. 

Free particle wave-functions. The free particle wave- 
functions in the angular momentum basis can be conve- 
niently expressed in terms of the following vectors: 

(z) = [Z m (z), Z m -i(z), Z m (z), Z m+1 (z)] , (7) 



w(p) 



\ 



+ V/2) 2 ~p 2 ](u;- 
[(u + V/2) 2 -p 2 ]p 
t ± (u 2 - V 2 /A) 
t±(u-V/2)p 



V/2)\ 



/ 



(8) 



The last vector is useful as long as ui 7^ V/2. The func- 
tion determining the eigenstates is: D(p,u>) = [p 2 — 

V74 - u 2 ] 2 + V 2 t\/A - u 2 (V 2 + t\). Then, provided 
that D(k,u>) = 0, (k > 0) which corresponds to prop- 
agating modes, the eigenfunctions are proportional to 
*z,m(w,fe,r) = Ui tm *Vz,m(kr)*w(k), where Z m (z) = 
J m (z) or Y m (z) are Bessel functions. The star product 
of two vectors is a vector with components defined by 
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[a ★ b] c 



If on the other hand D(iK,u)) 



0, 



(Re[«;] > 0) the eigenfunctions are: 



*jf,m(w,K,r) = Mo, m *%,mW *«)(£«), (9) 
*/ im (w,K,r) = 1i ,m *!)/,mH * w(-in), 

with I m (z) and K m (z) being modified Bessel functions. 
That these vectors are eigenstates can be verified directly 
by applying the real space version of |T|) to them. 

Local impurity wave-functions. The bound state wave- 
functions can be read off from the t-matrix equation. For 
a Dirac delta potential the resulting wave-functions are 
proportional to the propagator from the site of the im- 
purity to the site of interest evaluated at the energy of 
the bound state. These propagators are easily expressed 
in terms of the free-particle wavefunctions, e.g.: 



Gr a j,Al 
G aJ; Bl 



— i2M 2 A 2 



(10) 
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There are also analogous contributions coming from the 
other valley. The asymptotic behavior of the modified 
Bessel functions K n (z) ~ exp(— z)j^fz as z — > oo im- 
plies that the bound states are exponentially localized 
on a scale that is the same as in the general consid- 
eration above in ([5]) and J6]). At short distances one 
may use that K n (z) ~ l/z n for n > 1 to conclude 
that the wave-functions are not normalizable in the con- 
tinuum. The divergence is however not real since in a 
proper treatment of the short-distance physics it is cut- 
off by the lattice spacing a. The characteristic size of 
the wave-functions allows for a simple estimate of the 
critical density of impurities n c above which the overlap, 
and hence the interaction between the impurities, be- 
comes important. For weakly bound states we estimate 
n c ~ {Vtj_/(k g t)] 2 /(2TrV3) {E h /E g ), indicating that the 
critical density increases with the applied gate voltage. 
For V < £_l we find n c » 2.5 x 1(T 3 (E h /E & ), and for 
U < 1 eV one has n c < 10~ 6 . This result shows that even 
a small amount of impurities can have strong effects in 
the electronic properties of the BGB. 

Variational approach. A simple variational wave- 
function consisting of a wave-packet with angular mo- 
mentum m and a momentum close to k a 



from the E^ 



band is: * var (£) cx fg* dp*j tm (E + -(p),p,r) s /p/t, 
where we assume that £ -C k g . A variational calcula- 
tion shows that for any m, a weak attractive potential of 
strength oc U leads to a weakly bound state with binding 
energy E^ oc U 2 . This can be understood by noting that 
for each value of m, due to k g being non-zero, the problem 
maps into a ID system with an effective local potential, 
and in ID a weak attractive potential (cx U) always leads 
to a bound state with binding energy Eh cx U 2 . Thus the 
result is a direct consequence of the peculiar topology of 
the BGB band edge 



Potential well. The free-particle solutions can be used 
to study a simple BGB potential well modeled by two 
potentials: gj = —jjQ(R — r)/R. 9 is the Heaviside step 
function so that R is the size of the well and jj its di- 
mensionless strength in plane j. Bound states of the po- 
tential well are then described by two ^x,m outside and 
the appropriate pair of \&j, m 's and ^j^'s inside of the 
well. By matching the wave-functions at r = R we have 
studied the binding energies and find that the deepest 
bound states are in one of the angular momentum chan- 
nels m = 0,±1 for a substantial parameter range. Since 
these types of states are also present for the Dirac delta 
potential we argue that the physics of short-range poten- 
tials can be approximated (except for the short-distance 
physics) by Dirac delta potentials with a strength tuned 
to give the correct binding energy. A typical result for 
the binding energies is shown in Fig. [TJ The important 
case of a screened Coulomb potential generally requires 
a more sophisticated approach. Nevertheless, we do not 
anticipate any qualitative discrepancies between a poten- 
tial well and a screened Coulomb potential. We expect 
the screening wave-vector to be roughly proportional to 
the density of states at the Fermi energy; and once the 
range and the strength of the potential have been esti- 
mated a potential well can be used to estimate the bind- 
ing energies. We also note that the asymptotic behavior 
in ([5]) is quite general for a decaying potential. 

Coherent potential approximation. As discussed above, 
for a finite density m of impurities the bound states inter- 
act with each other leading to the possibility of band gap 
renormalization and the formation of impurity bands. A 
simple theory of these effects is the CPA [22|, [23| . In this 
approximation the disorder is treated as a self-consistent 
medium with recovered translational invariance. The 
medium is described by a set of four local self-energies 
which are allowed to take on different values on all of the 
inequivalent lattice sites. The self-energies are chosen 
so that on average there is no scattering in the effective 
medium. Explicitly we introduce the diagonal matrix 
H-e(oj) = Diag[SAi, Sbi, ^A2, £b2] , here and in the fol- 
lowing we suppress the frequency-dependence of the self- 
energies for brevity. Then the Green's function matrix is 
given by 



G-^w.k) 



Ho(k)-W s (w). 



(11) 



Finally, following the standard approach to derive the 
CPA jlil . [23|], we obtain the self-consistent equations: 
T, a j = mU/\l — (U — E Q j)G a jl. Using the propagators 
obtained from (fTT|) it is straightforward to compute G us- 
ing the same approximations that lead to j3|) . From these 
equations one can obtain the density of states (DOS) on 
the different sublattices: p a j(u>) — —Im.G a j(u> + i5)/Tt. 
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In the clean case one finds: 
- V/2 



P°Al 



P°B1 



2A 2 



X ■ 



uV(2 - x) 



^{V 2 + t 2 ± )u 2 -VH'i/A 
t\(uj + V/2)(2-x) 



(12) 



for \uj\ > EJ2. Here X = (0,1,2) for (|cj| < V/2, 

V/2 < \u\ < V*± + ^ 2 /4, y/*± + V2 / 4 ^ M)- The 
corresponding quantities in plane 2 are obtained by the 
substitution V — > — V. Takingthe limit V — > we re- 
cover the unbiased result of Ref.Lfl. Notice the square-root 
singularity that starts to appear already above V/2. 
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FIG. 2: (Color online) Left: DOS as a function of the energy 
(in units of E g ) close to the conduction band edge for differ- 
ent impurity concentrations (see inset), U = —1 eV. Right: 
Details of the DOS inside of the gap for different impurity 
concentrations for U — > oo. In both cases V = 40 meV. 

The numerically calculated density of states for U — > 
oo is shown in Fig. O The impurity band evolves from 
the single-impurity B2 bound state which for the param- 
eters involved is located at e w 0.3E g . Further evidence 
for this interpretation is that the total integrated DOS 
inside the split-off bands for the two lowest impurity con- 
centrations is equal to rij. It is worth mentioning that 
the width of the impurity band in the CPA is likely to 
be overestimated. The reason for this is that that the 
use of effective atoms, all of which have some impurity 
character, facilitates the interaction between the impuri- 
ties (23| . For smaller values of the impurity strength the 
single-impurity bound states are all weakly bound (cf. 
Fig. [1]) and the "impurity bands" merge with the bulk 
bands as shown in Fig. [2] and [3] The bands have been 
shifted rigidly by the amount UiU for a more transparent 
comparison between the different cases. The smoothen- 
ing of the singularity as well as the band gap renormaliza- 
tion is apparent. Observe also that the band edge moves 
further into the gap at the side where the bound states 
are located. It is likely that the CPA gives a better ap- 
proximation for these states since by ([5]) they are weakly 
damped almost propagating modes. Notice that the gap 
and the whole structure of the DOS in the region of the 
gap is changing with V, and in particular the possibil- 
ity that the actual gap closes before V = because of 
impurity-induced states inside of the gap. 




FIG. 3: (Color online) DOS as a function of energy (in units 
of Eg) for different values of the applied bias V (see inset) 
and U = -10 eV. Left: n; = 10" 3 ; Right: m = 10 -2 . 



Conclusion. We have studied the effects of impurities 
in a biased graphene bilayer. We find that local poten- 
tials always generate bound states inside of the gap. A 
finite density of impurities can lead to the formation of 
impurity bands as well as band gap renormalization. We 
show that, unlike ordinary semiconductors, the electronic 
properties and density of states can be controlled exter- 
nally via an applied voltage bias. 

We thank A. Geim, F. Guinea and N. M. R. Peres for 
many illuminating discussions. A. H.C.N, is supported 
through NSF grant DMR-0343790. 
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